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Abstract  —  This  paper  presents  theoretieal  results  on 
instability  proeesses  that  oeeur  in  the  positive-differential- 
resistanee  region  of  nanoseale  tunneling  structures  and 
reports  on  efforts  to  development  advanced  numerical 
techniques  for  use  in  future  optimization  studies.  These 
results  were  obtained  from  numerical  implementations  of  the 
Wigner-Poisson  electron  transport  model.  Here,  the  primary 
focus  of  the  reported  research  is  on  developing  simulation 
methods  that  are  adaptive  to  parallel-computing  platforms. 
Together,  these  investigations  demonstrate  the  high 
computational  demands  associated  with  modeling  fully  time- 
dependent  phenomenon  in  resonant  tunneling  structures 
(RTS)  and  offer  new  numerical  solutions  for  the  rapid  and 
efficient  analysis  of  these  types  of  problems.  Furthermore, 
the  simulation  tools  under  development  will  enable  future 
investigations  into  new  quantum  phenomenon  that  strongly 
influence  instability  processes  in  RTSs. 


I.  Introduction 

Resonant  tunneling  diodes  possess  extremely  fast 
response  times  and  have  been  utilized  as  the  gain  element 
in  oscillator  sources  at  frequencies  approaching  1 
terahertz  (THz).  Unfortunately,  the  traditional 
implementation  of  an  RTD  oscillator  (i.e.,  where  the  RTD 
is  used  as  an  extrinsic  gain  element)  leads  to  low  output 
powers  that  are  on  the  order  of  microwatts  [1].  The 
limited  output  power  of  the  extrinsic  RTD  oscillator  is 
directly  a  result  of  the  broadband  gain.  Specifically,  when 
an  RTD  device  is  embedded  in  an  electrical  circuit  or 
cavity,  there  will  be  a  tendency  to  induce  current 
oscillations  across  the  entire  frequency-band  where  the 
gain  exists.  Hence,  high-frequency  output  power  will 
either  be  lost  to  unwanted  spurious  low-frequency  modes 
or  significantly  scaled  back  due  to  reductions  in  device 
area  (and  therefore  output  current)  that  must  be  used  to 
lower  RTD  capacitance  to  achieve  low-frequency 
stabilization.  These  basic  facts  have  motivated  a 
theoretical  search  for  an  intrinsic  RTD  oscillator  design 


that  utilizes  a  microscopic  instability  mechanism  directly 
[2].  In  fact,  our  most  recent  efforts  have  been  towards  the 
specification  of  a  resonant  tunneling  structure  (RTS)  that 
will  yield  current  oscillations  within  the  positive- 
differential-resistance  region  of  the  average  current- 
voltage  (I-V)  characteristic.  This  particular  scheme  holds 
promise  because  it  has  the  potential  to  realize  a  device 
with  an  intrinsic  instability  mechanism  that  will  not 
produce  oscillation  modes  outside  a  narrow  prescribed 
domain.  Hence,  it  may  be  possible  to  practically 
implement  such  an  intrinsic  RTS  oscillator  circuit  without 
reducing  the  device  contact  area  (along  with  the  output 
current  and  power)  as  was  needed  for  the  low-frequency 
stabilization  of  the  extrinsic  RTD  oscillator. 

Very  recent  simulation  studies  performed  upon  emitter 
engineered  RTSs  suggest  that  it  may  be  possible  to  induce 
current  oscillations  within  the  PDR  region  of  the  average 
I-V  characteristic  [3].  Specifically,  a  Wigner-Poisson 
model  was  utilized  to  study  the  electron  dynamics  of  an 
RTS  with  an  emitter  engineered  region  as  shown  in  Fig.  1. 
Here,  the  fractional  compositions  of  the  Gai.yAlyAs 


Fig.  1.  Details  of  an  emitter-engineered  (EE)  RTS. 
material  in  the  emitter  region  have  been  graded  to 
prematurely  induce  an  emitter  quantum-well  (EQW) 
within  the  undoped  portion  of  device  in  front  of  the  first 
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barrier.  This  was  done  to  alter  the  multi-subband  strueture 
that  eouples  the  EQW  to  the  main  quantum-well  (MQW) 
of  the  RTD  under  a  eondition  of  applied  bias.  The  goal  of 
this  partieular  engineering  is  to  eneourage  a  eertain  type 
of  subband-eoupling  between  the  EQW  and  MQW  whieh 
have  been  shown  previously  [3]  to  be  the  underlying 
instability  meehanism  for  intrinsie  oseillations  in  RTDs. 
The  resulting  average  I-V  eurves  for  the  RTS  strueture 
defined  in  Fig.  1  are  given  in  Fig.  2.  Note  the  emergenee 


Fig.  2.  Average  I-V  results  for  the  EE-RTS  in  Fig.  1. 

of  two  bias  voltage  windows  (BVWs)  with  one  residing  in 
the  PDR  region.  In  addition,  there  is  an  enhaneement  in 
the  amplitude  of  the  eurrent-density  oseillation  in  the 
seeond  BVW  as  shown  in  the  results  given  in  Fig.  3. 


eireuit.  However,  these  simulations  of  the  isolated  RTS 
deviee  are  very  eomputationally  intensive  and  required 
over  a  week  of  runtime  on  a  serial  platform.  Obviously, 
the  full  time-dependent  simulation  of  the  RTS  deviee 
embedded  in  an  oseillator  eireuit  will  require  the  use  of 
effieient  numerieal  tools  on  parallel  eomputing  platform. 
This  paper  will  report  on  efforts  to  develop  advaneed 
numerieal  teehniques  for  use  in  future  optimization 
studies.  Here,  the  emphasis  is  on  the  development  of  fast 
and  effieient  tools,  that  are  adaptive  to  parallel  eomputing 
platforms,  and  that  ean  be  used  for  future  investigations 
into  new  quantum  phenomenon  that  influenee  the 
operation  of  RTS-based  oseillator  eireuits. 

I.  WiGNER-POISSON  ELECTRON  TRANSPORT  MODEL 

The  Wigner  funetion  formulation  of  quantum 
meehanies  was  seleeted  for  these  investigations  into  the 
RTS  beeause  of  its  many  useful  eharaeteristies  for  the 
simulation  of  quantum-effeet  eleetronie  deviees,  ineluding 
the  natural  ability  to  handle  dissipate  and  open-boundary 
systems.  The  Wigner  funetion  ean  be  eombined  with  the 
Poisson  equation  to  provide  for  an  adequate  quantum 
meehanieal  deseription  of  the  eleetron  transport  through 
tunneling  nanostruetures.  Details  regarding  the  derivation 
ean  be  found  elsewhere  [4],  but  the  model  is  a  two 
equation  system  with  the  basie  mathematieal  form 
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Fig.  3.  Oseillation  results  for  the  seeond  BVW  in  Fig.  2. 


preliminary  results  suggest  that  it  may  be  possible  to 
design  RTS-based  oseillator  eireuits  that  do  not  induee 
extraneous  low-frequeney  oseillations  in  the  external 
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where  the  last  term  in  Eq.  (1)  is  given  by 
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with  the  physieal  eonstant  Cj  =  —  hj(2zm*)  and  the 
integral  expression  defined  by 
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where  L  is  the  length  of  the  tunneling  structure  under 
consideration.  The  last  term  in  Eq.  (3)  is  due  to  scattering 
dissipation  and  is  modeled  using  the  relaxation  time 
approximation  [1].  The  boundary  conditions  on 
f  {x, k' ) at  the  emitter  ( x  =  0  )  and  collector  (x—L)  are 
specified  to  approximate  flat-band  transport  and 
equilibrium  electron-distribution  conditions  [5]. 

III.  Temporal  Integration  Algorithms 

A.  Integration  Method  -  To  track  the  time-evolution  of  the 
current  density  for  a  given  value  of  applied  bias  voltage, 
the  domain  and  equations  are  discretized  and  given  to  a 
numerical  integrator,  which  accurately  solves  for  the 
Wigner  function,  f,  for  given  time  interval  over  the 
discretized  domain.  So  for  a  given  time  t,  this  method 
computes  a  vector / such  that  (fly  ~  f(Xi,kj,t),  which  is  the 
descretized  Wigner  function  over  space  and  momentum 
coordinates.  These  approximations  are  computed  over  the 
coordinates,  x,-,  i^l,2,...,N^  and  kj,  j^l,2,...,Nh  Two 
separate  methods  were  investigated  and  compared  in 
terms  of  efficiency.  While  designing  appropriate  doping 
profiles,  barrier  well  placement  and  widths,  applied 
voltage  biases,  and  other  controllable  physical  parameters, 
several  temporal  simulations  will  be  required.  If  the 
method  is  too  inefficient,  some  design  methods  could  be 
impractical  to  perform  or  even  infeasible.  The  first 
integrator  considered  was  ROCK4  [6].  ROCK4  is  a  fourth 
order  Runge-Kutta  method  which  incorporates  variable 
stages  to  enhance  the  method's  stability.  ROCK4  is  an 
explicit  method,  meaning  the  current  solution  calculated 
by  the  integrator  is  based  only  on  previously  computed 
data.  Therefore,  the  advantage  of  this  method  is  that  it 
requires  low  storage  and  the  solution  can  be  quickly 
calculated  from  known  values.  The  disadvantages  of 
these  methods  are  their  stability  regions.  In  order  for  the 
algorithm  to  remain  stable,  an  explicit  method  is  forced  to 
take  small  time  steps,  forcing  the  integrator  to  do  more 
work.  The  second  method  explored  by  the  authors  was 
VODEPK  [7,8].  This  method  is  a  variable-order  implicit 
method.  Since  it  is  implicit,  at  each  time  step  a  nonlinear 
equation  must  be  solved  to  compute  the  current  solution. 
While  this  is  a  disadvantage  when  compared  to  explicit 
methods,  the  implicit  methods  are  allowed  to  take  larger 
time  steps  since  these  methods  are  not  concerned  with 
satisfying  stability  constraints.  So  if  the  nonlinear 
equations  could  be  solved  effectively,  then  the  implicit 
method  could  have  an  overall  advantage  over  the  explicit 
method. 


B.  Preconditioner  Development  -  To  solve  the  nonlinear 
equations  that  arise  in  the  time  integration,  VODEPK  uses 
a  Newton-GMRES  algorithm.  This  is  an  iterative  method 
that  finds  the  root  of  a  nonlinear  function 

ffl  ffl 

G  :  R  — >  R  where  the  next  iterate  is  found  by 
stepping  forward  from  the  current  iterate/,.  Note,  for  our 
problem,  and  f„  is  the  initial  guess  for  the 

Wigner  function  at  the  current  time.  The  step  direction, 
s„,  solves  the  linear  equation  G’(f„)s„^-G(fiJ,  where  G’(f„) 
is  the  Jacobian  matrix  of  G  evaluated  at/,.  Therefore,  a 
linear  equation  must  be  solved  for  every  Newton  iterate. 
Here,  the  equation  is  solved  with  an  iterative  method 
called  GMRES.  This  method  does  not  require  the 
computation  of  the  coefficient  matrix  (i.e.,  Jacobian 
matrix  in  our  case)  or  its  storage,  which  is  advantageous 
since  as  we  refine  the  grids  to  get  a  more  accurate 
simulation  of  the  RTD,  the  dimension  of  this  problem  is 
expected  to  become  quite  large.  Also,  it  has  been  well 
established  that  if  a  good  preconditioner  can  be  found, 
GMRES  can  efficiently  find  solutions  to  large  systems  of 
linear  equations  [9].  A  preconditioner  is  another  matrix, 
M,  multiplied  into  the  linear  equation  in  hopes  of  better 
conditioning  the  Jacobian  matrix,  hence  accelerating  the 
convergence  of  GMRES.  So  GMRES  would  solve  the 
linear  equation  MG'(fn)s„^-G(fJ,  where  MG’(f„)  would  be 
a  better-  conditioned  coefficient  matrix  than  the  usual 
G'(fy).  The  Newton  iteration  successfully  terminates  when 
either  the  size  of  the  next  step  is  below  a  set  tolerance,  or 
norm  of  the  nonlinear  function  is  below  a  set  tolerance. 
Note  that  both  of  these  parameters  will  be  approaching  a 
very  small  value  if  the  algorithm  is  converging  to  the 
correct  answer.  When  applying  VODEPK  to  the  Wigner 
equation,  the  Jacobian  matrix  of  the  nonlinear  equations  in 
VODEPK  have  the  form  G'{f)  =  I  -  A  ■  dW{f)l  d  f  , 
where  A  is  some  constant.  If  we  ignore  all  the  terms  in 
W(f)  except  the  D^(f)  term,  then  we  get  that 
dW{f)ld  f  ~Dj^.  Therefore,  a  good  preconditioning 
matrix  would  be  M  =  (I  -  ADJ'^  since  G  « (7  -  ADJ'k 
This  preconditioner  is  also  computationally  inexpensive  to 
implement  since  the  matrix  representation  of  is  sparse. 

C.  Numerical  Results  and  Statistics  -  The  new  integration 
algorithms  were  implemented  and  tested  in  our  base 
simulation  code  [10]  used  to  derive  the  results  given  in 
Section  I.  It  is  important  to  note  that  the  original  code 
utilizes  a  coherence  length,  L„  that  is  tied  to  the  upper 
value  of  the  x-space  grid.  Hence,  it  is  also  important  to 
monitor  the  value  of  as  we  refine  the  grid.  The  table 
below  summarizes  the  serial  runtimes  of  the  ROCK4  and 
VODEPK  version  of  the  simulations  on  several  different 
grids.  Note  that  all  of  the  computations  were  performed  on 
the  IBM-SP3  at  the  North  Carolina  Supercomputing 


Center  (NCSC).  The  runtimes  below  suggest  that 
VODEPK  handles  time  integration  more  effieiently  than 
ROCK4,  espeeially  as  the  grids  are  refined.  So  VODEPK 
has  been  ehosen  for  our  future  simulation  sinee  its 
performanee  will  be  better  for  generating  aeeurate 
solutions  to  larger  dimensional  problems.  Figure  4  shows 
the  average  I-V  eurves  for  the  (N^,Nk)  grids  of  (86,72)  and 
(64,128).  Note,  the  (86,72)  grid  results  exaetly  duplieate 


Table  I.  CPU  Time  versus 


N. 

Nk 

ROCK4 

(min) 

VODEPK 

(min) 

Lc 

(A) 

Ax 

(A) 

Ak 

(m/A 

) 

32 

32 

10 

5 

550 

17.8 

5.5 

32 

64 

40 

60 

1118 

17.8 

2.8 

64 

64 

124 

108 

550 

8.7 

5.6 

86 

72 

294 

58 

459.4 

6.5 

6.7 

64 

128 

432 

365 

1109 

8.7 

2.8 

simulations  that  were  made  with  the  original  eode  -  see 
previously  published  results  in  [10].  The  (86,72)  grid 
results  in  Fig.  4  are  in  almost  exaet  agreement  with  those 
generated  earlier  with  the  original  eode;  however,  the  new 
simulator  obtained  the  results  approximately  100  times 
faster.  Note  also  that  the  results  in  Fig.  5  show  that  the 
overall  results  are  strongly  influeneed  by  the  value  of  Lc, 
Flenee,  future  work  will  foeus  on  investigating  the  effeets 
that  both  the  eoherenee  length  and  the  applied  bias  voltage 
has  on  ereating  and  sustaining  eurrent  oseillations  in 
engineered  RTSs.  Furthermore,  the  physieally  appropriate 
value  of  Lc  will  be  determined  through  eomparisons  to 
experimental  results. 

IV.  Conclusion 

The  results  of  these  eomputational  studies  have 
demonstrated  an  advaneed  numerieal  algorithm  for  the 
fast  and  effieient  study  of  resonant  tunneling  struetures. 
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